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ABSTRACT 


'A finite element or Galerkin type semidiscrete method is proposed for numerical 
solution of a linear hyperbolic partial differential equation. The question of 
stability is reduced to the stability of a system of ordinary differential 
equations for which Dahlquist theory applies. 

We also present some results of separating the part of numerical solution which 
causes the spurious oscillation near shock-like response of semidiscrete scheme 
to a step function initial condition. In general all methods produce such 
oscillatory overshoots on either side of shocks. This overshoot pathoTogy , 
which displays a behaviour similar to Gibb's phenomena of Fourier series, is 
explained on the basis of dispersion of separated Fourier components which 
relies on linearized theory to be satisfactory. We present expository results, 
polished formal proofs will appear elsewhere. 


INTRODUCTION 


Our model of one and two dimensional linear hyperbolic equations are 
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Introducing bS C -cCcsot and Cu-cSia<X equation ( 2) can be written as 

X J 


(2') '^ + cc rt oi , aJ + cSu«|a=o 

2t 2>x ^ 

Galerkin or finite element semidiscretization [1], [2] , [3J , [4] .seeks an 
approximate solution for equation ( 2 * ) in the form 


( 3) 


^ t-t (t’) 

m.n. Jnva 0' Ml 

where 



( 4) 
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m n 10 Otherwise. 


We obtain a system of ordinary differential equations by requiring that the 

residual te orthogonal to the basis functions '■p^i.e 

7x. 'ly 

^ t R.^»=0* Candidates for are too many producing algorithms with 
increasing complexity proportional with their smoothness. We only present bilinear 


finite elements on squares. The orthogonality requirement yields, say in one 
dimensional case 
(5) 

where K^and are discrete Toeplltz operators with eigenvectors n V . 

if K w is an identity operator then scheme is explicit, otherwise implicit. If 
the real part of the corresponding eigenvalue is zero then the scheme is 
conservative [6] , [7] . The quantity 

(6) Z (jx>) - - Jn i 

u) 

Is the velocity of propagation of numerical solutions in comparison with exact 
propagation velocity C In ( 1). The quotient C(iu)/cor difference Co**)- C in 
an appropriate norm is the measure of spurious oscillations and dispersions In 
numerical solutions. Purely mathematical treatment without the effects of 
discretization i.e nonnumerical can be found in [8j. 

In the next sections, to study the response of semidiscrete scheme to sharp 
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gradient changes we simulate a shock by a step function Initial condition In ( 1), 

here we present an heuristic argument for the cause of parasitic oscillations 

around a point of discontinuity. 

Consider the weighted Galerkin semidiscretization 

(7) s£ C_ („ _ u \ 

2 at 2. -3t ■ 2H V 

of our model equation ( 1), where ot6lo,tl is a parameter. Note that ct*o corresponds 
to the equation 

( 8) centered difference approximation to 

* 

Since for any n , In equation ( 7), indices take three successive Integer 
values we may relabel them for n even as u.^ and for n odd as we then 

obtain respectively the following systems 

( 9) 
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dt*' = ” ^ 
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These equations are consistent approximations for the following systems 
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V . i 

A 


t d« 


Eliminating u or V In (11) we obtain respectively 


( 13 ) 
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Showing that in a doubly spaced grid wave equation Is satisfied. This Indicates 
that finite differencing is consistent with (13) rather than ( 1). Also adding 
the equations in (11) 



we see that discretization Is consistent for the avarage of the solutions at 
two successive grid points. However subtracting equations in (11) we obtain 
(15) 

This shows that due to discretization difference, however small, of two successive 
solutions propagates as an error wave in the discrete medium in the opposite 
direction. 

For (12), adding we obtain 



and subtracting we find 

in) 

which is the cause of oscillations In general. 


GALERKIN SEMIDISCRETIZATION 

FOR EQUATION 

(2’) 

On the square with vertices 

^ ^nt-1 > y r\+i 

>- (’WV.I K 1 >ni 

< x *., . v. 1 

we take basis 

functions to be 

1 

Vv 

for 

^ ~ ^ <- x - 


\ - 

K 

for 


(18) <fU, 3 )= < 
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for 

x ' X rn-,| 


Vv 

for 



0 otherwise 
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These are pyramids whose base- is a square with vertices are given above, centered 
at (X ,u ) with unit height. Forming the Inner products with the residual we 
obtain 

(19) v - fc - 

Only nonvanishing terms come for the values of indices ks m-t^m^rvw-iand L« * 

T t hus equation (19) reduces to 
\ \ 

(20) .IP . . r u ,^-w^a^c.a . B = O 

fer, ^ ‘Xwk,<a — 3 ^- -^ 7 r J — + d m-k/v-S. — 7 

Computation of inner products as double integrals are straightforward but 
tedious. Replacing the values of various integrals in equation (20), we obtain 

{21 } srtw-l + 4 U A.,M+ V"-' + ,*n“\ 

where fl= c *lS ot and ^=^£2 . 

6 J r c, ' *• 3 ^ 3 

system of equations (21) can be written in matrix notation on a rectangle [ 0,(M+l)h^ 

X [0,(N+1 )h in various ways. Let . . u^] T , (is 4,1,.. H , I N be the 


where 


N*N identity matrix, and L N =[ ^Cjl Nx ^ 

i J ' r l -< 

l i ] 

LO otherwise, 

superscript T indicates transposition. Then 

H « + 41 n + L ^) u v + ( L ->‘ + 4I « + L - > = _ kW-' L » +L « j + 

(.r L N +ci xlN+ e O'J 2 ]-i fi ^+^ + ^-(eL f j+^ x i N 4-(iL T N )vj 0 ^- 

BC + AI h +U n )U 0 + 


* (22) ^^l( L N +4I N +L I) tJ «-> + 4(L.,+41 N +L) l )0^+(.L N -AI+\- K )'J ni( J»: 
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*i^t' , ;rv-,- i - w ^+ ?< -»v + vl-^ for k.<». 


■MVJ 


3G>it *-Ii)UaV 

~ 2Vul“ ( e Ur'"V^-M + ^ ^~h) ^m-|" V ^)t) M V k,^M* W *A**» 

where we Introduced vectors In IR,^ 

<% 

v l = t- 6 V'°-- '°'-rMw1 r . fc«°A- •••*-» 

V >°>- v 0 ■ , fc«l,2,...,»A 

Xk ’tf tU fc.>°'---' °- eU fc,M+\] T , k = 2,3,...,KW\ 

^fc.“L u, ito> °> ' ' ' -* 0 > “•' ,x+i 1 , k 50 ,l,...|Mt|, 


We let 
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Then the system (22) in vector and block tridiagonal matrix notation becomes 


±.L 

36 dt 
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where 


*T 

• ■ ’3«-ra^v>i 


Here entries of matrices are NXN matrices and entries of vectors are N vectors. 
Note that for time Independent boundary conditions the last term in this equation 
vanishes. Further simplification is obtained by introducing NM dimensional vectors 
or M dimensional compound vectors i.e vectors whose components are N dimensional 
vectors, 
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and the square matrices of order NM, \A for the matrix on the left and (£> for 
the matrix on the right hand side of equation (23). 

The linear system (22) or equivalently (23) can be written as 

(24) SMJU-C 
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with the initial condition 

when t= 

0 

This system has a unique solution for 4U 

A, = 4 

r T <* 

\ 0 

T 

0 
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•n 
T«0. 


direct multiplication shows that and A 2 commute, this Is a direct con- 
sequence of both being Tceplltz matrices, therefore they have the same eigenvectors. 
Eigenvalues of A^, as easily verified, are 


nM-v-1 

with corresponding eigenvectors 






nh+\ 


Nt*W\ J 


Let be an eigenvalue of A2 associated with the eigenvector since A[ 

we have 

ftX x, ,X fc > “^PvXfc'Xj^ s < A ^w*X> ~ 1 XlX? * pvcX^»X^> 

and X^O Implies so . Gerschgorln theorem applied to A 2 yields 

ft. + V >0 ’ hence ^ ls nonsin 9ular [9 }. It Is known that A Is unltarlly 
similar to a dlagoial matrix £) with eigenvalues of A , which are the sum of the 
eigenvalues of A^ and A 2 , are the dlagona 1 entries. This similarity transfor- 
mation Is performed by taking Ss^X. > . . . .^'X Jl l.e columns of S are 

-i -\ 

eigenvectors of A . Letting S \JUU and multiplying (24) by S the Initial 
value problem reduces to 

B^ttS''vASU]»^^-f)IU=-^. ii S' 1 SSIU +S''c 

SX«LI 

Note that one does not ne^d to compute S , since S = &*• 

For the solution of (25) one step methods such as Runge-Kutta method can be used. 

Also a large number of multistep methods, Implicit or explicit in time (predlctor- 

5-8 , 
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corrector), o.ce a starting procedure is realized by a one step method, are 
available, and their stability theory is well understood and detailed treatment 
can be found in [10] , [113 - 

To show that the finite differencing scheme is conservative, we must show that 
the eigenvalues of Galerkin difference operators in (21) belonging to 

eigenvectors exp 1 are purely imaginary where o,' x = u) and 

. Substituting = a (Jt) c ^ n 3 * n (21 ) after some 

manipulation yields for the left hand side 


+ e ik ^x-^ e 

J 


and for the right hand side 

+ ** ( e^i. 2 ”* ) ')> 
Lc x s «^ X W (^2 [X"*' CcowO^Vv^ 


Hence 

aja wo.x'i 

where 


bifttOyW \ 


C ■ ?•, \ 


(26) %(u >,*) = J oc I Cooot 5-1 r + ^ 

* L cu k 1 + lCxn^lx *V V l-CoSiO.W 


3 <7 3 


1 


Si J t 3 w “3' 

which is imaginary. Setting we find the numerical solution 

^•) = Oy C o) ex^> X [ oJ x x m + y a - u: c 0°;*H ] • 

The discrepancy between C(wi ( ot) and C or more precisely the order of zero of 
Cca>,«()-C about ujH*Q is the the order of accuracy of the semidiscrete method. 
To show that this method is of order four, we expand C^u^ct.) in a Taylor series 

1 and a straightforward computation shows that 

•| (27) C(u?,gl)-c =■ - (ujWf-v 6 (tu x wf 4 (y^-f ) • 
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ERROR ESTIMATES FOR PRE AND POST OSCILLATIONS ABOUT DISCONTINUITIES 


To justify the heuristic argument presented earlier we may assume that this 
spurious oscillations are represented by small perturbations in a> , and replace 
U) by In trial solutions 

(28) aCx,i)=a t e :iWr8Cx - ?w+ ^ :l 

Expanding C(u> r£) in a Taylor series about £-0 and retaining only the 
linear terms we obtain 

Since u>»£ , terms of order £. can be neglected and introducing group velocity 

( 29) §Cui) = (.uj coj^)) 

(28) can be written as 

(30) 

Straightforward computation shows that eigenvalues of ( 7 ) corresponding to 
eigenvectors ^€> l ' u;X ' v lj , are 

Su\.uiV\. 

l»oC-V-oLCoiu 3k 

ar.d therefore 

C(<a)- - - r— 

i-oL-t-ocC^iu; h 

Using (29), is easily computed as 

(31 ) Q (jji) = C 

d (j -tx-h <x CoSoiW.) x 

Due to the discretization of the domain of the equation, the group velocity 
corresponding to 2 h wavelenght, from equation (31) is 

(32) 

which is the same as depicted in (17) and for oL^c in (15). 

To obtain estimates on local and global error of numerical solution we recall 
the definitions. [5] p.43, [13 J. We say an infinite series T 4 is (C,l) 




V 




.V" 


^k*iL u. » in this case 
tuo '• 


summable if lim (X -Jim ^±^£L^irL^ = s exists wher^ 
we write Z^* 6 ^ c » 1 ) sense. 

An infinite series Z 1 ^ is said tobe summable by Abel's method ( some say 
Poisson s) or simply A-summable to s, if is convergent for \r\<\ and 

rSiX 14 -!^ Where S fe is defined above. 

We need two results, the first is that the series 


00 


(3J) — - -4- Cosvvx 
^ n=\ 

is (C,l) and also A-summable to zero. 

It, Is known that, [5J p.20 
S ■= *4- ■+■ Coskx = — (,A V «| x 

fcrl blA.=. 

and from the trigonometric identity 
2-5U| Ccskx- CGSCt-HW- 

S; f ? ^ 


3^-i 


4ii 


,rr 


it follows that 

Si = 

fc-O 

Thus 

n. 


SuV 


C. * WrtA-y 

0«.a x - 


~ U-V\ 21 ~ rv-v \ 2. Sin.:* 21 ^ Ln -C^ + iV x - 

. ^ Kro 2L 

n-Jc?** 0 • 10 show A-summability, recall the Poisson's formula [5] , p.61; 

l+ir n CW= 


2L 


^2-lrr 


\-r 

ZO-ZrCo^x-f-r 1 ) 


lr\<\ ' 


Letting r-*l we see that the assertion is true. 


The second result is that 

tO 


(34) "*} Suvrtx 

is (C ,1 ) and also A-summable to _L ^ oc 

2. 2. 

It is known that [5], p.21; 

T VX * 

S s Y Si<\(vx c: — Cot — — t Ssitklhj L aV x ' 

ETi X 2. i £; (V k 

Using the trigonometric identity 

Z SlaJ. \Co<> ^ ’V- Cos • • • +G.s(a'vi)xU 6in.C v '--v'Vx - Sln.x 
we find 
Or - 

»V “ * * 

h.4-1 


. r 


i 


= -ii— -i Cot i 

1 1 ASJi 


s 


and the result follows by letting n->oo 

5-11 


S'«g(rit \)?Zz. SGvx. 
a V \ 




To show A-summability we use Poisson's formula 

Z ji c. r &uv.x 

r J1.1 ax =■ 


n- 

and let r 


A-2rCa5x.-+ r- 


\. 


He now estimate the X^norm of the global error. As a direct consequence of 

Parseval's identity, it is known that Fourier transform is an isometric isomorphism 

between the Hilbert spaces Involved [16] p.51-52, [15] p.25. Therefore It 

suffices to compute the i^norm of the Fourier transform of the error. To simulate 

the shock, we let the initial condition to be the step function 
\ x£o 
O oc<c 

Without loss of generality we may assume that the discrete Fourier transform of 


(35) Ui.x 




the net initial condition u^(.o) is equal to the Fourier transform of (35) , and 




we obtain 

a. A -uox.i i v -ux>n.Vv. 

UrO » 

It follows from the proofs of statements concerning equations (33) and (34) that 


CO 


(36) UC^O>«a tv (oJ J o , ) = k21u.^o-)e 


series on the right of equation (36) is (C,l) and hence A-summable to 

a >lgj> 

(37) U«.w,0)=U.^U,o)= 

Zc bwi -j- 

From equation (1), Fourier transform of the exact solution is easily computed 

A a -icuifc 

U(u>;t) = Ut^ y o) e 

U a (£) is the solution of semidiscrete equation, for simplicity we assume 
to be the identity operator, taking the discrete Fourier transform of the 
semidiscrete equation and solving the resulting differential equation one 


obtains 

A 




For conservative schemes >ioj)=-uucW) , therefore the norm of the global 


error is !L 
K 


Z . f’ v A 2 -U*>CW)t -Lojct l 

jn- J \^>°)\ |e -e. \ 


^ir 

K 


^':V: 


\>’v ' ^ 




i \ 

;v. 


. I Introducing dimensionless variables "c. = t£ and u=u)K, > a straightforward 

k J 

• computation shows that 



» 

' CONCLUSION 


The semidiscrete method proposed here has a reasonable Courant number and a fourth 

t 

order accuracy. Results are theoretically conclusive. Computational evidence for 
detailed comparison of this method with conventional methods will await our 
numerical experiments. 

The measure of oscillations in the numerical solution, in a neighborhood of sharp 
changes is the pointwise error. We were able to show with a lenghty argument, 
although there are some gaps in details of proofs, that maxima of the difference 
between the exact and the numerical solutions continually diminish and minima 
continually increase in an interval of lenght 4h on each side of the sharp 
gradient change. Numerical solution is approximately 0.28h larger in the 
upstream direction. 
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